#define fast_norm3( v ) (v).w = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z+(v).z); (v).x *= (v).w; (v).y *= (v).w; (v).z *= (v).w; (poot).w = 0.0f;
#define fast_norm3_t( temp, v ) (temp) = native_rsqrt((v).x*(v).x + (v).y*(v).y + (v).z*(v).z ); (v).x *= (temp); (v).y *= (temp); (v).z *= (temp);

#define fast_dot3( v1, v2 ) ((v1).x*(v2).x + (v1).y*(v2).y + (v1).z*(v2).z)
#define fast_sub3( r, v1, v2 ) (r).x = (v1).x - (v2).x; (r).y = (v1).y - (v2).y; (r).z = (v1).z - (v2).z;
#define fast_add3( r, v1, v2 ) (r).x = (v1).x + (v2).x; (r).y = (v1).y + (v2).y; (r).z = (v1).z + (v2).z;
#define fast_min3( r, v1, v2 ) (r).x = fmin((v1).x, (v2).x); (r).y = fmin((v1).y, (v2).y); (r).z = fmin((v1).z, (v2).z);
#define fast_max3( r, v1, v2 ) (r).x = fmax((v1).x, (v2).x); (r).y = fmax((v1).y, (v2).y); (r).z = fmax((v1).z, (v2).z);

typedef struct {
	float x;
	float y;
	float z;
} float3h;

#define INNER_LOOP_1 \
	normal_lookup			= read_imagef(normSpecImage, nearestSampler, coords);\
	average_normal.x		+= 	normal_lookup.x; \
	average_normal.y		+=  normal_lookup.y; \
	average_normal.z		+=	normal_lookup.z;\
	depth	= read_imagef(depthImage, nearestSampler, coords).x; \
	minDepth = fmin( minDepth, depth ); \
	maxDepth = fmax( maxDepth, depth ); \
	coords.y				+=  1.0f;\
		
#define INNER_LOOP_2 \
	normal_lookup	= read_imagef(normSpecImage, nearestSampler, coords); \
	normal.x		= normal_lookup.x - 0.5f; \
	normal.y		= normal_lookup.y - 0.5f; \
	normal.z		= normal_lookup.z - 0.5f; \
	fast_norm3_t(temp, normal); \
	min_dot = fmin( min_dot, fast_dot3(normal, average_normal));\
	coords.y		+= 1.0f; \

#define UNROLL2(x) x x
#define UNROLL4(x) UNROLL2(x) UNROLL2(x)
#define UNROLL8(x) UNROLL4(x) UNROLL4(x)
#define UNROLL16(x) UNROLL8(x) UNROLL8(x)
#define UNROLL32(x) UNROLL16(x) UNROLL16(x)

__kernel void perblockattributes(	__read_only	 image2d_t	normSpecImage,
									__read_only	 image2d_t	depthImage,
									const		sampler_t	nearestSampler,
									const		int2		imageSize,
									const		int2		num_blocks,
									__write_only image2d_t	perBlockImage1, 
									__write_only image2d_t	perBlockImage2)
{
	float temp;
	int i;

	int2 range;
	const int2 gid = { get_global_id(0), get_global_id(1) };
	range.x = min( BLOCK_X, imageSize.x - mul24(gid.x , BLOCK_X)); 
	range.y = min( BLOCK_Y, imageSize.y - mul24(gid.y , BLOCK_Y)); 

	if (range.x < BLOCK_X || range.y < BLOCK_Y) {
		return;
	}
	
	float2 baseCoords = { mul24(gid.x , BLOCK_X), mul24(gid.y , BLOCK_Y) };
	float2 coords;

	/*
		This code section computes the bounding box
	 */
	float maxDepth = 0.0f, minDepth = 1.0f;

	float normal_cone_size;
	float3h average_normal = { 0.0f, 0.0f, 0.0f };

	// compute average norm and min/max depth
	coords.x = baseCoords.x;
	for (i=0; i<range.x; i++) {
		
		float4 normal_lookup;
		float depth;
		
		coords.y = baseCoords.y;

#if BLOCK_Y == 1
		INNER_LOOP_1
#elif BLOCK_Y == 2
		UNROLL2(INNER_LOOP_1)
#elif BLOCK_Y == 4
		UNROLL4(INNER_LOOP_1)
#elif BLOCK_Y == 8
		UNROLL8(INNER_LOOP_1)
#elif BLOCK_Y == 16
		UNROLL16(INNER_LOOP_1)
#elif BLOCK_Y == 32
		UNROLL32(INNER_LOOP_1)
#else
	#error INVALID BLOCK_Y VALUE
#endif

		coords.x += 1.0f;

	}
	average_normal.x -= ( BLOCK_X * BLOCK_Y ) * 0.5f;
	average_normal.y -= ( BLOCK_X * BLOCK_Y ) * 0.5f;
	average_normal.z -= ( BLOCK_X * BLOCK_Y ) * 0.5f;
	fast_norm3_t(temp, average_normal);
	// end compute average normal
	
	// find cone size
	float min_dot = 1.0f;
	
	coords.x = baseCoords.x;
	for (i=0; i<range.x; i++) {
		
		coords.y = baseCoords.y;
		
		float4 normal_lookup;
		float3h normal;
		
#if BLOCK_Y == 1
		INNER_LOOP_2
#elif BLOCK_Y == 2
		UNROLL2(INNER_LOOP_2)
#elif BLOCK_Y == 4
		UNROLL4(INNER_LOOP_2)
#elif BLOCK_Y == 8
		UNROLL8(INNER_LOOP_2)
#elif BLOCK_Y == 16
		UNROLL16(INNER_LOOP_2)
#elif BLOCK_Y == 32
		UNROLL32(INNER_LOOP_2)
#else
	#error INVALID BLOCK_Y VALUE
#endif
		
		coords.x += 1.0f;
		
	}
	// get normal cone radius squared from minimum dot product
	normal_cone_size = min_dot <= 0.0f ? -1.0f : 1.0f - min_dot * min_dot;

	/* compute ndc versions of min and max depths */
	float maxZ = ( 2.0f * maxDepth - 1.0f );
	float minZ = ( 2.0f * minDepth - 1.0f );	
		
	float4 out1 = { average_normal.x, average_normal.y, average_normal.z, normal_cone_size };
	float4 out2 = { minZ, maxZ, 0.0f, 0.0f };
	
	write_imagef( perBlockImage2, gid, out2 );
	write_imagef( perBlockImage1, gid, out1 );

}